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Abstract 



The energy spectrum of cosmic rays between 10 eV and 10 eV, derived from 
measurements of the shower size (total number of charged particles) and the total 
muon number of extensive air showers by the KASCADE-Grande experiment, is 
described. The resulting all-particle energy spectrum exhibits strong hints for a 
hardening of the spectrum at approximately 2 • 10 16 eV and a significant steepening 
at ~ 8 • 10 16 eV. These observations challenge the view that the spectrum is a single 
power law between knee and ankle. Possible scenarios generating such features are 
discussed in terms of astrophysical processes that may explain the transition region 
from galactic to extragalactic origin of cosmic rays. 



1 Introduction 



The main goals of experimental cosmic ray research are the determination 
of the arrival direction distribution, the primary energy spectrum, and the 
elemental composition. Those measurements comprise important hints to un- 
derstand the origin, acceleration and propagation of energetic cosmic particles. 
The needed measurements can be done directly or indirectly, depending on the 
energy of the primary particle. At energies above 10 15 eV, the energy spectrum 
must be determined indirectly from the measured properties of extensive air 
showers (EAS) that cosmic rays induce in the Earth's atmosphere pQ. 

The determination of the primary energy and elemental composition in the 
energy range from 10 15 eV up to above 10 20 eV is subject of earth-bound ex- 
periments since more than five decades. It has been shown that the all-particle 
spectrum has a power-law like behavior (oc i?~ 7 , with 7 rs 2.7) with features, 
which are known as 'knee' and 'ankle' at 3-5-10 15 eV and 4-10-10 18 eV, respec- 
tively. Whereas at the knee the spectrum steepens, the ankle is characterized 
by a flattening of the spectrum by roughly the same change of the spectral 
index of A7 = ±0.3-0.4. Cosmic rays above the ankle are most probable of 
extragalactic origin j2], i.e. somewhere in the energy range from 10 16 eV to a 
few 10 18 eV a break-off of the heavy component and the transition of cosmic 
rays of galactic to extragalactic origin is expected. 

As the measured position of the knee is roughly in agreement with the energy 
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where supernova remnants (SNR) become inefficient accelerating particles [3], 
various theories with different assumptions were developed to explain the be- 
havior of the spectrum between the knee and ankle features. 
The basic idea of the 'dip model' [I] is that the ankle is a propagation feature of 
extragalactic protons (energy loss by electron pair production). Consequently, 
in that model the composition at the ankle is to a large extent proton-dominant 
and the transition from galactic to extragalactic origin of cosmic rays occurs 
already at energies well below 10 18 eV. In the scenario of the dip model, at 
energies around 10 17 eV a pure galactic iron component should be left and a 
small kink in the spectrum at around 5-7- 10 17 eV, as indicated by observations 
by the AKENO [5] and HiRes [6] experiments and named as 'second knee' [7], 
would be assigned to the transition. This is in agreement with the SNR the- 
ory, where the knee positions of individual primary masses are proportional 
to the charge of the nuclei starting with the proton knee at around -Ef nee = 3- 
5 • 10 15 eV and E^ nee = Z ■ E^ nee (rigidity dependence of knee positions for 
galactic cosmic rays). 

On the other hand, to avoid an early appearance of the extragalactic cosmic 
ray component, Hillas j3] proposed in addition to the standard SNR compo- 
nent, a 'component B' of cosmic rays of galactic origin. This component would 
also experience a charge dependence of break-offs, but now shifted to approx- 
imately ten times higher energy. As a result, the transition occurs here at the 
ankle and for the entire energy range from 10 15 eV to 10 18 eV a mixed ele- 
mental composition is expected. In this scenario, the second knee, if it exists, 
would be a feature of the component B. 

The KASCADE experiment and its extension, KASCADE-Grande, aim to 
provide high quality air-shower data in the energy range of 10 14 eV to 10 18 eV 
to evaluate the validity of these models and to distinguish between them. The 
KASCADE experiment has shown that the knee is due to a distinct break in 
the proton intensity despite protons are not the most abundant primary in 
this energy range. The break is followed by a kink in the spectrum of Helium 
nuclei [8], i.e. the knee in the all-particle spectrum is a feature of the light 
nuclei (Z< 6), only, where the difference in the energies of the knee features of 
primary protons and Helium facilitates the assumption of a charge dependence 
of the break-off. First analyses of KASCADE-Grande data jH] resulted in a 
knee-like feature at around 8 • 10 16 eV caused by a steepening in the spectrum 
of heavy primary cosmic rays. In the present analysis, the reconstruction of the 
all-particle energy spectrum of cosmic rays in the range from 10 16 to 10 18 eV 
is described in detail. 

Depending on the experimental apparatus and the detection technique of 
ground-based air-shower experiments, different sets of EAS observables are 
available to estimate the energy of the primary cosmic ray [10]. In case of 
ground arrays the total number of charged particles (often called shower size) 
in the shower and the corresponding particle density at observation level are 
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commonly employed. The muon content of EAS plays an important role, too. 
In the atmosphere the muon component suffers less attenuation than electro- 
magnetic or hadronic components and exhibits less fluctuations compared to 
the more abundant electromagnetic component. In KAS CADE- Grande both 
components, the muon and the electromagnetic ones, are measured with in- 
dependently operating detectors. Both, together with the information of their 
correlation on a single-event-basis, are used to derive the spectrum. After a 
short description of the apparatus and the reconstruction procedures of the 
EAS parameters, we will describe the method developed to determine the 
all-particle energy spectrum including studies of systematic uncertainties. We 
conclude this paper with a discussion of the results. 



2 The Experiment 

The experimental layout, as well as the reconstruction procedures and accura- 
cies of KASCADE-Grande observables are described in detail in reference [TTj . 
In this chapter, we only summarize the most important facts relevant for the 
present analysis. 

2.1 KASCADE-Grande 

The multi-detector experiment KASCADE [TJ] (located at 49.1°N, 8.4°E, 
llOma.s.l.) was extended to KASCADE-Grande in 2003 by installing a large 
array of 37 stations consisting of 10 m 2 scintillation detectors each (Fig. [TJ. 
KASCADE-Grande provides a sensitive area of about 0.5 km 2 and operates 
jointly with the existing KASCADE detectors. Main parts of the experiment 
used for the present analysis are the Grande array spread over an area of 
700 x 700 m 2 , and the original KASCADE array covering 200 x 200 m 2 . The 
Grande array is installed over an irregular triangular grid with an average 
spacing of 137 m. The KASCADE array is composed of 252 detector stations 
on a square grid with 13 m spacing. It is organized in 12 outer clusters of 16 
stations each and 4 inner clusters of 15 stations each. The outer clusters (192 
stations) are equipped with two unshielded (e.m.) and one shielded (muon) 
detector units each. A muon detector unit consists of 4 plastic scintillators of 
90 x 90 x 3 cm 3 each, where the iron-lead shielding provides a threshold of 230 
MeV kinetic energy for vertically incident muons. The total sensitive area of 
the muon array amounts to 622 m 2 . 

While with the Grande array we reconstruct the total number of charged 
particles, data from the shielded scintillation detectors of the KASCADE ar- 
ray are used to reconstruct the total number of muons on an event-by-event 
basis for events triggered by Grande. The combination of both allows us to 
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Fig. 1. Layout of the KASCADE-Grande experiment: The KASCADE array and 
the distribution of the 37 stations of the Grande array are shown. The outer 12 
clusters of the KASCADE array consist of shielded ^-detectors (hatched area). The 
dotted line shows the fiducial area selected for the present analysis. 

reconstruct the energy spectrum of cosmic rays in the range from 10 16 eV up 
to 10 18 eV. 



2.2 Monte Carlo simulations 



The simulations performed include the detailed air shower development in the 
atmosphere, as well as the response of the detector and its electronics. There- 
fore, the EAS parameters reconstructed from simulated showers are obtained 
exactly in the same way as for real data. The EAS were generated with uni- 
formly distributed core positions (at an area larger than the Grande array), 
with isotropically distributed arrival directions, and with a spectral index of 
7 = —2, i.e. roughly one order of magnitude harder than the measured spec- 
trum. The spectral slope is chosen as a compromise between computing time 
and statistical accuracy at the highest energies. Later, in the analysis proce- 
dures the simulated data is weighted to describe a softer energy spectrum with 
7 = -3. The EAS were simulated with CORSIKA [13] and the Monte Carlo 
generators FLUKA [2] and QGSJet II [15] (hadronic interactions). Sets of 
simulated events were produced in the energy range from 10 15 eV to 10 18 eV 
for five different representative mass groups: H, He, C, Si and Fe with about 
353, 000 events per primary. A few showers for all primaries were also sim- 
ulated for the higher energy range up to 3 ■ 10 18 eV in order to study the 
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reconstruction quality even beyond the energy range of interest. In addition, 
with less statistics, similar simulation sets were generated based on the high- 
energy hadronic interaction model EPOS, version 1.99 [TB] . 

2.3 Reconstruction 

Basic shower observables like the core position and the angle-of-incidence, as 
well as the total number of charged particles are provided by the reconstruction 
of data recorded by the Grande stations. The angle-of-incidence is determined 
using the relative particle arrival times at the stations. The core location, 
the slope of the lateral distribution function and the shower size (i.e. the to- 
tal number of charged particles N c h) are calculated by means of a maximum 
likelihood procedure, comparing the measured number of particles with the 
one expected from a NKG-like lateral distribution function [17] of charged 
particles in the EAS. KASCADE-Grande provides the unique opportunity of 
evaluating the reconstruction accuracies of the Grande array by a direct com- 
parison with an independent experiment. For a subsample of events collected 
by the Grande array it is possible to compare on an event-by-event basis 
the two independent reconstructions of KASCADE and Grande. By means of 
such a comparison the Grande reconstruction accuracies of the total number 
of charged particles are found to be: systematic uncertainty < 5%, statisti- 
cal accuracy better than 15%; accuracy of the arrival direction of about 0.8°; 
accuracy of the core position about 6 m. The total number of muons is 
calculated using the core position determined by the Grande array and the 
muon densities (E^ > 230 MeV) measured by the KASCADE muon array de- 
tectors, is derived from a maximum likelihood estimation comparing the 
measured densities with a lateral distribution function, where the uncertainty 
is less than 12%. For the purpose of the present analysis, the estimated 
has been corrected for systematic uncertainties using a correction function, 
which is derived from Monte Carlo simulations as explained in appendix IA.1I 
The correction improves the accuracy to better than 10%, including the un- 
certainty due to hadronic interaction models and unknown composition. 

2.4 Data selection and shower size spectra 

The selection used for the present analysis requires that the events passed suc- 
cessfully the full KASCADE-Grande reconstruction procedure. In addition, we 
requested events to be from stable periods of data taking with more than 35 
Grande stations and all 16 KASCADE clusters in operation. The selected 
events had to pass a 7/7 Grande hardware trigger (six of a hexagonal shape 
and the central one) and had to hit more than 11 Grande stations in total. 
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Fig. 2. Differential shower size N c h and muon number spectra for different zenith 
angular ranges. 



We restricted ourselves to events with zenith angles less than 40° due to in- 
creasing shower size uncertainties. Only air showers with cores located in a 
central area of KASCADE-Grande (see Fig. [1]) were selected. With this cut 
on the fiducial area, border effects are discarded and the under- and over- 
estimations on the muon number for events close to and far away from the 
center of the KASCADE array are reduced. All of these cuts were applied also 
to the Monte Carlo simulated events to study their effects on the selection, 
to optimize the cuts, and to control the uncertainty of the acceptance. Full 
efficiency for triggering and reconstructing air-showers, as well as a uniform 
distribution of the shower cores in the fiducial area, is reached at a primary 
energy of about 10 16 eV, slightly depending on the zenith angle, and on the 
primary particle type. 

The analysis presented here is based on 1173 days of data. The cuts on the 
fiducial area and zenith angle result in a total acceptance of 1.98 ■ 10 5 m 2 -sr, 
and an exposure of 2.0-10 13 m 2 -s-sr, respectively. Approximately 1.5-10 6 events 
are subject of the following analysis. 
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Fig. [2] shows the reconstructed N c h and distributions for five different 
angular ranges. The angular ranges were chosen in order to have the same 
acceptance in each of them. 



3 The reconstruction of the all-particle energy spectrum 



KASCADE-Grande, with the possibility of an independent reconstruction of 
the two observables N c h and per individual EAS, allows us to go for a 
dedicated strategy in estimating the all-particle energy spectrum. By means 
of Monte Carlo simulations, calibration formulas are obtained to calculate the 
primary energy per individual event based on N c h and N^, taking into account 
the correlation between the two particle components. This reduces the compo- 
sition dependence of the energy assignment. To account for attenuation effects 
in the atmosphere, which are different for the two shower observables, the en- 
ergy calibration is performed separately for the five zenith angular ranges. 
Finally, the obtained energy spectra are unfolded to account for the bin-to- 
bin migrations before they are combined to the resulting all-particle energy 
spectrum. 
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Fig. 3. Left panel: Scatter plot of the reconstructed N^/N^ vs. N c h for primary iron 
and proton nuclei and for the first angular bin. The full dots and error bars indicate 
the mean and statistical errors of the distribution of the individual events (small 
dots). The fits result in parameters c and d of expression [3l Right panel: Scatter 
plots of E vs. N c h for iron and proton primary nuclei. The fits result in parameters 
a and b of expression [TJ 
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3.1 Energy calibration 



The energy assignment starts by applying E = f(N c h, k), where k is defined 
through the ratio of N c /j and N^: k = g(N c h,N^). The main aim of the k 
parameter is to correlate these observables on an event-by-event basis, taking 
into account the differences in the N^/N^ ratio for different primary masses 
with the same N^, as well as the shower-to-shower fluctuations for events of 
the same primary mass: 



log w (E/GeV) = [a H + (a Fe - a H ) ■ k] ■ log w (N ch ) + 
+ b H + (b Fe -b H )-k 



k 



logioiNch/Nj - log w (N ch /N^ H 



log^Nch/N^pe - logioiNch/N^u 

log^Nch/N^H^e = C H ,Fe ' log W (N ch ) + d H ,Fe- 



(1) 

(2) 
(3) 



By definition the k parameter is a number centered around zero for proton ini- 
tiated showers and around one for iron initiated showers |18j . The coefficients 
a, b, c, d are obtained independently by simulations for each zenith angular 
range and for each primary mass, where fits are applied to the scatter plots 
(NchjNch/N^) and (N ch ,E). The fit range is chosen to be 6 < \og(N ch ) < 8, 
i.e. where 100% trigger efficiency is guaranteed. Primary protons exhibit larger 
fluctuations than heavier primaries, therefore, the coefficients c and d are ob- 
tained iteratively in case of protons in order to improve the reconstruction of 
the energy spectrum. As an example, Fig. [31 shows the scatter plots including 
the resulting functions for the first angular bin. Shown are the errors on the 
mean, which are small due to the large Monte Carlo statistics. For the fits, 
however, we also take into account the width of the distributions in order 
to avoid a bias due to varying shower-to-shower fluctuations, in particular in 
case of primary protons and small shower sizes. It is obvious that taking into 
account the correlation of the observables reduces significantly the composi- 
tion dependence of the energy assignment. Similar procedures are applied to 
the other angular bins and all the coefficients are compiled in Table [TJ The 
uncertainties of these numbers are small, but considered in the calculation of 
the total systematic uncertainty. 



Applying the derived energy calibration to the measured data, all-particle 
energy spectra for the five zenith angular ranges are obtained (Fig. |4]). To 
refine the energy assignment function from the so far assumed pure power-law 
behavior of the {N c h, N c h/ Ny) and (N c h,E) relations to a more realistic non- 
linear calibration, as well as to unfold bin-to-bin migrations due to shower- 
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Table 1 

Coefficients of the energy calibration functions. 



Angular bin 


a 


b 


c 


d 




H 


Fe 


H 


Fe 


H 


Fe 


H 


Fe 


< 16.7° 


0.91 


0.88 


1.33 


1.82 


0.10 


0.16 


0.79 


-0.06 


16.7° < 6 < 24.0° 


0.89 


0.88 


1.50 


1.92 


0.08 


0.18 


0.88 


-0.25 


24.0° < 6 < 29.9° 


0.94 


0.89 


1.30 


1.94 


0.10 


0.16 


0.68 


-0.17 


29.9° < 6 < 35.1° 


0.93 


0.88 


1.46 


2.10 


0.11 


0.17 


0.54 


-0.35 


35.1° < 6 < 40.0° 


0.92 


0.88 


1.75 


2.29 


0.11 


0.16 


0.41 


-0.35 



to-shower fluctuations, response matrices Rij for the different angular bins 
are constructed and applied, i.e. the spectra are unfolded (see appendix I A. 2|) . 
Effects of this procedure (see Fig. HJ) on the flux are estimated to be smaller 
than 5% at all energy bins and, therefore, do not significantly change the shape 
of the spectra. For the following discussions we always refer to the unfolded 
spectra. 

The spectra of the different angular ranges exhibit small systematic shifts 
relative to each other (see Fig. HJ), where we observe a slight flux increase with 
increasing zenith angle. This corresponds to a horizontal shift in the energy 
assignment, which could be explained if real showers penetrate deeper in the 
atmosphere than predicted by the QGSJetll hadronic interaction modePl. 
As the effect is observed for different assumed composition models and for 
all hadronic interaction models, it is most probably caused by a mismatch 
between the predicted and measured attenuation lengths of the shower particle 
numbers N ch and (for a more detailed discussion see appendix IA.3j) . The 
aforementioned differences in the spectra are considered as one of the major 
sources of the systematic uncertainty on the energy spectrum, and are taken 
into account in the estimation of the total systematic uncertainties. The final 
all-particle spectrum of KASCADE-Grande is obtained (see Figs. El d El and 
Table [3]) by combining the spectra for the individual angular ranges. Only 
those events are taken into account, for which the reconstructed energy is 
above the energy threshold for the angular bin of interest (see Fig. H]) . 

3.2 Systematic uncertainties 

Different sources of systematic uncertainties, which affect the all-particle en- 
ergy spectrum, are investigated. Most of the effects lead to a shift in energy 
and this, in turn, to a shift in the spectrum (see also Table [2]): 

• Attenuation: The average difference between the intensities obtained for the 

5 In fact, the analysis of muon production heights with the muon tracking detector 
of KASCADE leads to a similar conclusion, see ref. |19j . 
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Fig. 4. Reconstructed all-particle energy spectra for all five angular bins. In the 
upper panel the direct reconstructed as well as the unfolded spectra are displayed, 
where the spectra are scaled for better visibility. In the lower panel only the unfolded 
spectra are shown without scaling. Only statistical uncertainties are displayed. 

various angular bins has been used to define the systematic uncertainty as- 
sociated with the angular dependence of the parameters appearing in the 
energy calibration functions of the different angular ranges. This includes 
the systematic uncertainty related to the description of the air shower atten- 
uation in the atmosphere for simulated data. The uncertainty is used as well 
at lower energies, where not all angular bins contribute to the spectrum. In 
addition, this uncertainty indicates the limit for which the QGSjetll model 
reproduces the shower attenuation in a consistent way for the selected data 
sample. 

• Energy calibration and composition: To estimate a possible bias in repro- 
ducing the energy spectrum, eqns. [1] and [2] have been applied to simulated 
energy spectra build up by pure H and Fe primaries as well as the other 
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Table 2 

Estimated uncertainties of the cosmic ray intensity for different energies, where only 
absolute values are given in case of symmetric uncertainties. 



Source of uncertainty 


10 16 eV (%) 


10 17 eV (%) 


10 18 eV (%) 


intensity in different angular bins (attenuation) 


-0/+6.5 


10.9 


21.3 


energy calibration and composition 


10.3 


5.8 


13.4 


slope of the primary spectrum 


a n 


9 n 


1 Q 


reconstruction (core and shower sizes) 


0.1 


1.4 


6.5 


total 


-11.1/+12.8 


12.6 


26.1 


artificial spectrum structures (extreme cases) 




<10 




hadronic interaction model (EPOS-QGSJet) 


-5.3 


-16.9 


-14.6 


statistical error 


0.6 


2.7 


17.0 


energy resolution (mixed composition) 


24.7 


18.6 


13.6 
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Fig. 5. Left panel: Ratio between the reconstructed and initially simulated energy 
spectrum for light (blue), heavy (red) and mixed primaries (black) summing up all 
angular bins. The ratio using simulated EPOS data when the spectrum is recon- 
structed using the QGSJet-II based calibration functions is also shown. The larger 
fluctuations are due to the lower statistics of the available EPOS simulations. Right 
panel: Resolution of the energy assignment for simulated sets (QGSJet) with mixed 
composition, pure proton, and pure iron primaries, respectively. The dots show the 
offsets of the reconstructed energy in bins of simulated energy for mixed composi- 
tion, while lines represent the results for pure H and Fe assumptions. The circles 
and dotted lines show the corresponding RMS of these distributions. 

three mass groups (He, C, Si), and by a uniform mixture of the five pri- 
maries with 20% abundance each. The spectra have been composed to fall 
off with slope 7 = -3 and are compared with the reconstructed ones to find 
the corresponding uncertainty in intensity. Figure [5] summarizes the results 
on the ratio between the reconstructed and the initaliy simulated spectra for 
the full angular range. The original energy spectra are fairly well reproduced 
in all energy bins within a systematic uncertainty smaller than 10%. This 
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is quite important because it guarantees that any primary spectrum with 
arbitrary composition can be fairly well reproduced, even though the param- 
eters of the response matrix and of the muon number correction functions 
have been defined on the basis of a mixed composition. Another impor- 
tant aspect to investigate is the capability and robustness of the analysis to 
recognize changes of the spectral shape and the elemental composition in 
the considered spectrum. Different cases have been simulated, from rigidity- 
dependent knees with variable abundances for each single chemical group to 
completely artificial compositions, such as sudden changes around 10 17 eV 
from nearly pure iron to nearly pure protons and vice versa. Even for these 
extreme cases, the ratio between the reconstructed and initially simulated 
intensity deviates from unity throughout by less than 10%. 
Spectral slope of the Monte Carlo simulations: A further source of system- 
atic uncertainty is the choice of a spectral slope of 7 = -3 in the simulations 
to determine the energy calibration functions and response matrices. There- 
fore, new calibration functions and response matrices have been calculated 
using simulated spectra with 71 = -2.8 and 72 = -3.2. The difference be- 
tween the intensities obtained with the new coefficients has been defined 
as systematic uncertainty and results into 4% at low energies, where the 
fluctuations are more important, and decreases to about 2% at Ef^lO 17 eV. 
Reconstruction quality of shower sizes: The Grande array has an asymmet- 
ric geometry. Therefore, we checked the reconstruction of the observables 
by determining the energy spectrum for different ranges of distance of the 
shower core to the muon detector. The relative difference in intensity as a 
function of energy is used to compute a systematic uncertainty induced by 
the reconstruction, and it amounts to about 2% at around 10 17 eV, slightly 
increasing with energy. 

Energy resolution: Simulated data using an equal mixture of all primaries 
have been divided in bins of true energy and the distributions of the relative 
differences between reconstructed and true energies have been compared. 
As shown in Fig. El right panel, the RMS of such distributions (energy 
resolution) is about 25% at lower energies and decreases to about 15% at the 
highest energies due to the decrease of intrinsic shower fluctuations. Results 
for pure H and Fe primaries are also indicated by lines. As expected, proton 
initiated showers show larger fluctuations compared to EAS generated by 
primary iron nuclei. 

Hadronic interaction models: By now, for all considerations the model com- 
bination QGSJet-II/FLUKA has been used. As the calibration depends on 
simulations, other interaction models may change the interpretation of the 
data. To study such effects we investigated the influence of the hadronic in- 
teraction model by performing the energy assignment based on simulations 
with the hadronic interaction model EPOS (version 1.99). Due to the smaller 
statistics of the EPOS simulations, larger uncertainties are obtained and no 
response matrix corrections could be applied. Therefore, a proper energy 
spectrum by means of EPOS simulations cannot be derived yet, but general 
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characteristics are inferred. Comparing the all-particle energy spectra for 
both cases, it was found that EPOS leads to a slightly lower intensity (~10- 
15%) compared to QGSJet. A similar conclusion can be drawn, if the data 
simulated with EPOS are treated like experimental data and reconstructed 
using the QGSJet based calibration functions (see Fig. [5]). Now, as a con- 
sequence, the intensity is reconstructed 10-15% higher than the simulated 
input. The main reason behind the difference between QGSJet and EPOS 
results has to be ascribed to the different N^/N^ ratio predicted by the two 
models. In particular, EPOS predicts that the showers are richer in muons 
and slightly poorer in charged particles (for the relevant energy range, at 
sea level, and for the experimental conditions of KASCADE-Grande). These 
results are used as a rough estimate of the uncertainties due to the choice 
of the particular interaction model QGSJet-II. 



3.3 Cross-checks of the energy spectrum reconstruction 

Using the N^-N^ ratio we reduced the dependence of the reconstructed all- 
particle spectrum on the elemental composition. But, since both observables 
are reconstructed independently, we can apply an energy spectrum reconstruc- 
tion on both observables individually. At the end one would expect the same 
result for the energy spectrum by all approaches, provided that (i) the mea- 
surements are accurate enough, (ii) the reconstructions work without failures, 
and (iii) the Monte Carlo simulations describe correctly the shower develop- 
ment and its fluctuations, and (iv) the composition is known. But, the fact that 
not all of the above requirements are fulfilled and the individual observables 
exhibit substantial differences in their composition sensitivity hampers such 
straightforward cross-checks. Nevertheless, such analyses are used to check the 
reconstruction procedures and the influence of systematic uncertainties. Some 
details of these analyses [36|37] and the results can be found in appendix IA.4I 

In figure |6]the spectra obtained by the three methods are compiled, where the 
intensity is multiplied by a factor of E 3 ' . Owing to the different approaches, 
the results using single observables are given only for the pure proton and iron 
assumptions, whereas the final spectrum is displayed with a band showing the 
systematic uncertainties. Using N c h as the observable obtained from data of 
the Grande array only, a large dependence on the primary elemental compo- 
sition is present, which is reflected in a big difference between the intensities 
for proton and iron assumptions. The muon number shows less compo- 
sition dependence compared to the shower size, though it is still the largest 
contribution of uncertainty. The narrower range for a solution provided by 
compared to N c h confirms the finding of KASCADE that at sea-level the 
number of mostly low-energy muons A^ is less composition sensitive than the 
total number of charged particles [2D]. The method finally applied for energy 
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primary energy [eV] 

Fig. 6. Reconstructed all-particle energy spectrum by the three different approaches 
applied to KASCADE-Grande data. For the combined method a band is shown, 
indicating the range of the systematic uncertainty (without uncertainties due to the 
chosen hadronic interaction model). 

estimation, owing to the combination of the two variables, results in a larger 
reconstruction uncertainty. But the total uncertainty (including composition 
dependence) is considerably smaller by taking into account also the correla- 
tion of these observables. This additional information is strikingly decreasing 
the composition dependence. 

Of particular interest is the fact that by using N c h, the iron assumption re- 
sults in a higher intensity than the proton assumption, whereas using the 
opposite is the case (see appendix IA.4I) . In any case, if there is only the possi- 
bility of applying a one-dimensional method, then one finds a large variance in 
possible solutions (any solution within the range spanned by the proton and 
iron line, not even parallel to these lines). Interestingly, over the whole energy 
range there is only little room for a solution satisfying both ranges, spanned 
by N c h and A^, and this solution prefers a more heavy composition - in the 
framework of the QGSJet-II hadronic interaction model. The results obtained 
by combining N c h and A^ lie within the area spanned by the other methods. 
This finding expresses a cross-check of the intrinsic consistency of the results 
for the interpretation of two measured observables based on the used hadronic 
interaction model QGSJet-II. This was found to be valid also for the case of 
interpreting the data with the hadronic interaction model EPOS-1.99. 



4 Results and discussion 

The resulting all-particle spectrum exhibits structures which do not allow us 
to describe the spectrum with a single power law. To emphasize this, figure [7] 
shows the residuals of the all-particle energy spectrum multiplied by a factor in 
such a way that the middle part of the spectrum becomes flat. The power law 
index of 7 = —2.92 ± 0.02 is obtained by fitting the range of \og 10 (E / eV) = 
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primary energy [eV] 

Fig. 7. The all-particle energy spectrum obtained with KASCADE-Grande. The 
residual intensity after multiplying the spectrum with a factor of E 2,918 and nor- 
malized with A is displayed as well as the band of systematic uncertainty. 



Table 3 

Differential intensity values of the all particle energy spectrum for the QGSJet-II 
based analysis. The first column of errors denotes the statistical uncertainty, the 
second column the systematic uncertainty. 



bin number Ener 


gy [cV] 


dl /dE± stat. ± syst. 
[m- 2 s- 1 sr~ 1 GeV- 1 




1 


1.11 


10 16 


(2.46±0.02±°;§ 7 ) ■ 


10~ 


15 


2 


1.41 


10 16 


(1.16±0.01±8;i|) • 


io- 


15 


3 


1.78 


10 16 


(5.49±0.03±°;«)- 


io- 


16 


4 


2.24 


10 16 


(2.76 ± 0.02±g;f |) . 


io- 


16 


5 


2.82 


10 16 


(1.40 ±0.01 ±0.08) 


■ 10" 


-16 


6 


3.55 


10 16 


(7.13 ±0.07 ±0.50) 


■ 10" 


-17 


7 


4.47 


10 16 


(3.69 ±0.04 ±0.38) 


■ 10" 


-17 


8 


5.62 


10 16 


(1.89 ±0.03 ±0.23) 


■ 10" 


-17 


9 


7.08 


10 16 


(9.52 ±0.17 ±1.24) 


■ 10" 


-18 


10 


8.91 


10 16 


(4.67 ±0.11 ±0.58) 


■ io- 


-18 


11 


1.12 


10 17 


(2.19 ±0.07 ±0.28) 


■ 10" 


-18 


12 


1.41 


10 17 


(1.05 ±0.04 ±0.16) 


■ io- 


-18 


13 


1.78 


10 17 


(5.41 ±0.26 ±1.10) 


■ io- 


-19 


11 


2.24 


10 17 


(2.64 ±0.16 ±0.61) 


■ 10" 


-19 


15 


2.82 


10 17 


(1.12 ±0.10 ±0.28) 


■ io- 


-19 


16 


3.55 


10 17 


(5.66 ±0.59 ±1.34) 


■ 10" 


-21) 


17 


5.01 


10 17 


(2.47 ±0.23 ±0.61) 


■ io- 


-20 


18 


7.94 


10 17 


(4.26 ±0.75 ±1.11) 


■ 10" 


-21 


19 


1.41 


10 18 


(3.90 ± 1.36 ±0.96) 


■ io- 


-22 



16.2 — 17.0. For the full energy range a statistical analysis reveals that the 
spectrum is not described by a single power law with a significance of 2.1a. 
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Just above 10 16 eV the spectrum exhibits a 'concave' behavior, which is sig- 
nificant with respect to the systematic and statistical uncertainties. This is 
true despite the fact that only vertical showers contribute to the spectrum in 
this energy range (see Fig. H]). This hardening of the spectrum is validated 
by several cross-checks, e.g., by efficiency correction of more inclined events 
based on simulations. A hardening of the spectrum is, on the one hand, ex- 
pected when a pure rigidity dependence of the galactic cosmic rays is assumed. 
Depending on the relative abundances of the different primaries one would ex- 
pect charge dependent steps (i.e. slope changes) in the all-particle spectrum. 
The gap in the knee positions of light primaries (proton, helium, and CNO 
group of Z = 1—8) and the heavy group can lead to a hardening of the 
spectrum [21]. On the other hand, there are also other possible astrophysical 
scenarios to get a concave behavior of the cosmic ray spectrum. In general, 
a transition from one source population to another one could also result in 
a hardening of the spectrum. In such the KASCADE-Grande result 

could be a first experimental hint to the 'component B' of galactic cosmic 
rays, as proposed by Hillas [3]. A possible scenario for the component B is 
discussed by Ptuskin et al. [22], where the maximum acceleration energy for 
different types of supernovae is considered, taking into account also their rel- 
ative abundances in our galaxy. This scenario can lead to an extension of the 
galactic component up to a few EeV. In addition, in this model the transition 
of CR origin from the standard type SN la to SN lib supernovae requires a 
hardening of the spectrum at 10 PeV. General galactic modulation can lead 
only to a very smooth change of the slope index over more than a decade in 
energy, but postulating a contribution of a nearby (single) source, sharpening 
the knee at a few PeV [23|2l] , would also require a hardening of the spectrum 
just above 10 PeV. It is interesting to note that recently the CREAM detec- 
tor (balloon experiment) has described such a hardening of the proton and 
helium spectra at much lower energies [25]; which by the authors is assigned 
to a possible change of the acceleration mechanism of cosmic rays. 

Another feature in the spectrum is a small break slightly below 10 17 eV. Ap- 
plying a second power law above 10 17 eV an index of 7 = —3.39 ± 0.07 is 
obtained. The indices of the two power-laws differ from each other by two 
standard deviations. Even taking into account extreme scenarios for the sys- 
tematic uncertainties, or applying more stringent procedures to calculate the 
significance an effect with > lcr remains. Fitting the spectrum with a function 
of two power laws intercepted by a smooth knee the energy of the break is 
assigned to log l0 (E/eV) = 16.92 ± 0.09, which is in nice agreement to the 
value obtained by analysing the raw-like (i.e. not corrected for reconstruction 
uncertainties) all-particle spectrum [5J. In [9] it was also seen that the break 
gets more significant when analysing a subsample of events where showers gen- 
erated by heavy primary particles are enhanced. The change in slope occurs 
at an energy where the charge dependent knee of the iron component would 
be expected (KASCADE QGSJet based analysis assigns the proton knee to 
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an energy of ~ 3 • 10 15 eV). The change of the spectral index of this knee-like 
feature is small compared to the first one, original well-known knee [8], what 
could be explained, when the iron component is not dominant around 10 17 eV. 
This again can happen in presence of a 'component B' of mixed composition, 
but a final conclusion is not possible without investigating the composition in 
detail. 

Both observed features were subject to detailed cross-checks. In particular, we 
investigated how far the applied unfolding procedure affects the spectrum. To 
build up the response matrix an energy spectrum and a particular composition 
has to be assumed. We investigated possible effects by assuming extreme cases 
and by using different unfolding methods. If one assumes a very abrupt change 
of the spectral slope and in composition for a given energy, the resolution of 
KAS CADE- Grande would indeed smear that out to a structure distributed 
over values of 0.3 — 0.5 in log 10 (E/eV) of the reconstructed spectrum, but 
still clearly visible. 

At higher energies the KAS CADE- Grande spectrum, in particular close to 
10 18 eV, where other experiments have claimed a 'second knee' [7J, suffers 
from missing statistics. 

Despite the fact, that the discussed spectrum is based on the specific hadronic 
interaction model QGSJet-II, there is confidence that the found structures of 
the energy spectrum remain stable. The analysis has shown that the applied 
procedure can reconstruct the total number of charged particles, as well as the 
total muon number sufficiently well, independently of the hadronic interaction 
model in use. But the energy calibration assumes that the QGSJet-II model 
provides the correct lateral distribution of the particles over the entire distance 
range (exceeding the geometrical size of KASCADE-Grande). First studies 
with an alternative method to reconstruct the energy spectrum via the particle 
density at a fixed distance give hints to systematic deviations [2(5] in the energy 
calibration of the observable. But, the spectral structures discussed above are 
also present in the results of these studies. 

Figure [S] compiles the KASCADE-Grande energy spectrum with results of 
other experiments. Despite the independent measurements and data analy- 
sis there is a good agreement with the results of the KASCADE experiment 
and others in the overlapping energy range at low energies. In particular, the 
concave behavior seems to be needed to connect the spectrum with the spec- 
tra obtained by other experiments at the knee region. At higher energies the 
KASCADE-Grande spectrum (QGSJet-II) results in a slightly lower inten- 
sity compared to earlier experiments, in particular GAMMA, AKENO and 
YAKUTSK. The strong peak-like structure below 10 17 eV as was claimed by 
the GAMMA experiment [27] is not confirmed by our results. At the highest 
energy accessible by the KASCADE-Grande experiment, where we suffer from 
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Fig. 8. Comparison of the all-particle energy spectrum obtained with KASCADE- 
Grande data based on the QGSJet-II model to results of other experiments. The 
band denotes the systematic uncertainties. An analysis based on EPOS 1.99 would 
result in a spectrum which is shifted downwards by approximately 10% in intensity. 

missing statistics, our result is in agreement with a single power law and with 
the spectrum reported by HiRes and, when taken into account also the sys- 
tematic uncertainties mentioned for the Auger result, with the Pierre Auger 
Observatory. 



5 Conclusion 



The main air-shower observables of KASCADE-Grande, shower size and total 
number of muons, are reconstructed with high precision and low systematic 
uncertainties. Applying various reconstruction methods to the KASCADE- 
Grande data the obtained all-particle energy spectra are compared as a way 
to cross-check the reconstruction, to study systematic uncertainties and to test 
the validity of the underlying hadronic interaction models. By combining both 
observables, the all-particle energy spectrum of cosmic rays is reconstructed 
in the energy range of 10 16 eV to 10 18 eV within an uncertainty in intensity of 
10-15%, based on the hadronic interaction model QGSJet-II. 

Correcting the spectra for reconstruction uncertainties and taking into ac- 
count the systematic uncertainties for all methods, the underlying hadronic 
interaction models (QGSJet-II/FLUKA) result in a consistent solution, inde- 
pendent on the observable used, i.e. the single shower sizes or the correlation 
between the different observables. Tests with the hadronic interaction model 
EPOS 1.99 have shown that there is a shift in the absolute energy scale when 
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interpreting the data with this model, but the shape of the spectrum with 
its structures stays preserved. Progress in improving the interaction models is 
expected in the near future by detailed analyses of the now available data of 
the Large Hadron Collider, LHC (see, e.g. 



The resulting spectrum is consistent, and in the overlapping energy range 
in a very good agreement, with results of the KASCADE, EAS-TOP, and 
other experiments (Fig. [8]). The all-particle energy spectrum in the range from 
10 16 eV to 10 18 eV is found to exhibit some smaller structures: In particular, 
a hardening of the spectrum is observed at 2 ■ 10 16 eV and a small break- 
off at around 8 • 10 16 eV. These features are used to discuss the astrophysics 
in the transition region from galactic to extragalactic origin of cosmic rays, 
where a final conclusion is not possible without detailed knowledge of the 
elemental composition in this energy range. However, amongst others, the 
model proposed by Hillas [3J, e.g., which assumes a second component of 
galactic cosmic rays in addition to the standard SNR component, can explain 
the observed features of the measured all-particle energy spectrum. 



A wealth of information on individual showers is available with KASCADE- 
Grande. This makes it possible to reconstruct the all-particle energy spectrum 
with high precision, as well as to investigate the elemental composition, to test 
hadronic interaction models, and to study cosmic ray anisotropies. All these 
studies are under way and further results are expected in the near future. 
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A Appendix 



A.l Reconstruction of the total muon number 



Due to the fact that the muon detectors are located at the fringe of the Grande 
array, the uncertainty of the reconstructed muon number grows is seen to in- 
crease with the distance of the KASCADE array to the shower core from 
about 5-10% at 250 m to 25% at 600 m for log 10 (A^) > 5.3. But, as the fea- 
tures of these inaccuracies are well understood, we correct the reconstructed 
muon number by a correction function calculated on the basis of Monte Carlo 
simulations. This function takes simultaneously into account the dependence 
of the Nfj, uncertainties on the zenith angle of the reconstructed air shower, 
on the distance of the core position from the KASCADE array in shower co- 
ordinates, and on the muon number. The correction function was found to be 
nearly independent from the composition of cosmic rays and from the hadronic 
interaction model used. The uncertainty after applying these corrections is less 
than 8% for total muon numbers above log 10 (A^) > 5.3 (Fig. IA.1| . i.e. in the 
range of full efficiency 
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Fig. A.l. Accuracy (AA^/A^ with AA^ 
number before (solid circles) and after (open circles) applying the muon correction 
function versus the simulated muon size (left panel) and the distance of the EAS 
core to the KASCADE center (right panel). The lines represent the mean values 
of the achieved A^ accuracy for different primaries after applying the correction 
function and the band, the corresponding one-sigma width for a mixed composition 
assumption. 
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A. 2 Response matrix and unfolding 



As the fluctuations in the energy determination are larger than the bin size 
of the aimed-for energy spectrum an unfolding procedure is applied. Using 
Monte Carlo simulations a response matrix is constructed for the energy in- 
terval log 10 (£'/GeV) = 6 — 9.5, i.e. covering the entire range where fluctuations 
can affect the energy spectrum. This matrix represents the conditional proba- 
bility, P(Ej\Ej rue ), of an event with true energy in the bin log 10 (-E* rne ) being 
reconstructed with energy log 10 (-&,). By means of the response matrix a sys- 
tem of simultaneous equations, n^ xp = J2iLi P(Ej\Ej rue )nl rue , is established 
between the distribution of measured events, n% xp , and the actual energy dis- 
tribution, n\ rue . The system is solved iteratively for the true distribution by 
means of the Gold-algorithm [29|30] . Statistical uncertainties on the data are 
taken into account during unfolding by multiplying the equation system by 
the error matrix C*j = 5ij / a{n e * p ) [3D]. The unfolding procedure involves also 
a van Cittert transformation of the equation system to guarantee the positive 
definiteness of the modified response matrix, which is a necessary requirement 
for the convergence of the method [29] . On the other hand, unfolded results 
are validated with another recursive method known as the Bayes-algorithm 
|31j . In general, both unfolding techniques perform in a stable way for the 
present application. 

To avoid the problem of having wild fluctuations when increasing the number 
of iterations in the procedure a regularization method is applied, consisting of 
smoothing the result of unfolding in a given step before using it in the next 
iteration [31] . Smoothing was also applied to the response matrix, to avoid the 
presence of artificial effects in the unfolded distribution, which could arise from 
random fluctuations in some entries of the response matrix as a consequence 
of the limited statistics of the Monte Carlo data sets. Different methods for 
smoothing are applied and compared in order to find the optimal parameters 
for the procedure. To smooth the spectrum best results are obtained with the 
353if<5-twice algorithni"^1|32j . For the response matrix, quadratic fits along 
the diagonals are performed in the region of full efficiency to interpolate data 
into the region of low statistics. It is worth to mention that among the several 
tests employed to verify the performance of the methods, it was checked that 
the algorithms do not produce artificial structures in the spectrum or hide 
peaks which could be significantly present in the data, besides it was searched 
for the agreement between the forward-folded and measured distributions. 

For a large number of iterations, both the Gold- and Bayes-algorithms con- 
verge to the same result in the energy interval of full efficiency and good statis- 



6 In each step of the iteration, k, the input spectrum, plotted as log 10 [n^ ] vs 
log 10 (£' J - rite ), is smoothed with the aforementioned algorithm. 
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tics. In addition, positiveness of the solution is observed. As one goes deeper 
in the number of steps, the weighted mean of the squared sum of statistical 
errors and systematic bias, i.e (1/iV) X)f (°~stat,i + a syst,i) I 'nf ue , decreases and 
becomes constant when the result converges. This parameter speaks about the 
quality of the final solution. Unfolding without smoothing was also tried, but, 
in general, poorer results were obtained. 

A. 3 Attenuation correction 

In order to cross-check the attenuation effects of the observables N c h and 
on the treatment in reconstructing the all-particle spectrum, we applied to 
the two estimated observable spectra a correction which is independent of 
Monte Carlo simulations. The Constant Intensity Cut Method (CIC method) 
is based on the assumption that the arrival direction distribution of cosmic 
rays is isotropic and that the cosmic ray intensity and composition changes 
monotonically with primary energy [33|34|35] . In this way, the intensity of pri- 
mary particles becomes a reference variable for the primary energy of cosmic 
rays independent of the zenith angle. To apply the CIC method the integral 
spectra, J(> N c h) and J(> iV M ), are calculated for all angular bins. Then, 
fixed frequency rates (integral intensities) are chosen in the range of maxi- 
mum efficiency and sufficient statistics. By this, attenuation curves for each 
intensity are built, where an interpolation between two adjacent points of the 
integral spectrum is applied. The evolution of the shower size (muon number) 
in the atmosphere is extracted from the attenuation curves. The data have 
shown that it is possible to use constant attenuation parameters for the entire 
energy range since the differences in the obtained parameters by fitting indi- 
vidual curves are smaller than the uncertainty. With the parameters obtained, 
the shower size (muon number) of an individual air shower can be corrected 
with N ch ^(9 re f) = N ch ^(9) exp [P{6 re f) — P(9)] to obtain the equivalent size 
at a given zenith angle of reference, re f- The reference angle is chosen to be 
the mean of the measured zenith angle distribution, which is found to be 20° 
and 22° for the shower size and muon number distributions, respectively. Due 
to the independent reconstruction of N c h and slightly different reconstruc- 
tion thresholds lead to different mean angles and hence, different reference 
angles [36|37j . Uncertainties due to intrinsic assumptions on energy indepen- 
dent shower-to-shower fluctuations, the assumed spectral index, and constant 
composition are estimated and are taken into account. In Fig. IA.2l the result- 
ing spectra as corrected by the CIC method are shown for the shower size and 
for the muon number in the range of full efficiency for the independent re- 
construction of the observables. Note that the bin sizes are selected according 
to the resolution of N c h and N^. Applying this procedure to simulated data 
(QGSJet-II) has shown that indeed the differences between simulation and 
data are small for showers below 40° and that the estimated uncertainty on 
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Fig. A. 2. CIC corrected size spectra with statistical error bars and number of events 
per bin. 



our all-particle spectrum is a conservative approximation. 



A. 4 Single parameter reconstruction of the energy spectrum 



To reconstruct the energy spectrum by the observables independently, we use 
the reconstructed shower size spectra given in Fig. [2] and apply first an attenu- 
ation correction to the observables. Then, the shower size per individual event 
is calibrated by Monte Carlo simulations under the assumption of correlations 
in the form E Q oc N^ h and E oc N^, respectively, and an assumed primary 
composition [36f37] . 



To determine the calibration function of the number of (attenuation corrected) 
charged particles N c h and primary energy, Monte Carlo simulations were used, 
where a zenith angle range of 17° < 6 < 24°, i.e. around the reference angle, 
was selected. Assuming a linear dependence in logarithmic scale: log 10 E = 
a + b ■ log 10 N c h, the correlation between the primary energy and the number 
of charged particles is obtained, where the fit is applied in the range of full 
trigger and reconstruction efficiencies (see left panel of Fig. |A73|) . The fit yields 



a = 1.23 and b = 0.93 for primary protons and a = 1.75 and b = 0.90 for iron 
primaries. The numbers confirm the assumption of a power law and show that 
the slopes of these power laws are very similar. However, the normalization 
parameters assigning the energy are very different for the two primaries. The 
fits are also performed for helium, carbon, silicon, and for a uniform mixed 
composition to examine the dependence of the calibration on the assumed 
primary particle type, where the values of the fit parameters are found to 
be in between the above values. The energy resolution is estimated from the 
difference between true and reconstructed energy, resulting in 32% and 18% 
at E = 10 17 eV for proton and iron primaries, respectively, with an energy 
dependence of approximately l/y/E. 

To determine the energy from the (attenuation corrected) total muon number 
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Fig. A. 3. Calibration functions for assumed pure proton and iron primaries for the 
observables N c h (left panel) and (right panel), respectively. 
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Fig. A. 4. All-particle energy spectrum (including systematic uncertainties) for the 
assumption of primary protons and iron, respectively, based on the shower size (left 
panel) and on the total muon number (right panel). 

A^ a similar procedure as in the case of N c h was applied: a calibration function 
for the muon size (9 re f = 22°) in terms of the primary energy was invoked. As 
before, the calibration curve is described with a relation of the form log 10 E = 
a [l + b [l -\og w (see Fig. lA.3[ right panel). For the case of pure protons, the fit 
results in = 1.61 and = 1.09 and for iron primaries the values = 1.63 
and = 1.07 are obtained. Again it is observed that the slopes are similar, 
but here the difference in the parameter a is much smaller than in case of 
the charged particles. In case of the muon shower size, as the fluctuations are 
smaller, the energy resolution for protons and iron nuclei are of the order of 
25% and 12% (E w 10 17 eV), respectively. 



Figure IA.4I shows the all-particle energy spectra as obtained after applying 
the calibration functions, as well as the corrections for the bin-to-bin fluctua- 
tions for both approaches. Apart from the statistical uncertainties also bands 
are shown describing the systematic uncertainties on the intensity, which were 
determined for the two approaches and different primaries independently. The 
considered sources of such uncertainties include the estimate of the calibration 
functions, the chosen reference angle for the calibration, the muon number cor- 
rection function, and the application of the unfolding procedures. The total 
systematic uncertainty (i.e. sum in quadrature of all terms) on the intensity 
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for proton and iron is 21% (7%) and 10% (13%), respectively in case of N c h 
(iVjJ, at energies of 10 17 eV. The systematic uncertainties are energy depen- 
dent and evolve in such a way that they slightly increase near the threshold, 
where fluctuations are larger, and in the high-energy region, where statistics 
decreases. 
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